function gdr=Reg(Theta,Ni,regscheme)

%Computes the gradient of the regularization term of the wavelet transform
%of a velocity field Theta. The Ni's are matrices of the differential
%operators in the wavelet domain, regscheme is the scheme for
%regularization.

switch regscheme
    case 'HS'
        gdr=-cat(3,Ni(:,:,2)*Theta(:,:,1)*Ni(:,:,1)'+Ni(:,:,1)*Theta(:,...
            :,1)*Ni(:,:,2)',Ni(:,:,2)*Theta(:,:,2)*Ni(:,:,1)'+Ni(:,:,1)...
            *Theta(:,:,2)*Ni(:,:,2)');
    case '2Tik'
        gdr=cat(3,Ni(:,:,2)*Theta(:,:,1)*Ni(:,:,1)'+Ni(:,:,1)*Theta(:,...
            :,1)*Ni(:,:,2)',Ni(:,:,2)*Theta(:,:,2)*Ni(:,:,1)'+Ni(:,:,1)...
            *Theta(:,:,2)*Ni(:,:,2)');
    case 'lap'
        gdr=cat(3,Ni(:,:,3)*Theta(:,:,1)*Ni(:,:,1)'+Ni(:,:,1)*Theta(:,...
            :,1)*Ni(:,:,3)'+2*Ni(:,:,2)*Theta(:,:,1)*Ni(:,:,2)',Ni(:,:,...
            3)*Theta(:,:,2)*Ni(:,:,1)'+Ni(:,:,1)*Theta(:,:,2)*Ni(:,:,...
            3)'+2*Ni(:,:,2)*Theta(:,:,2)*Ni(:,:,2)');
    case 'divcurl'
        gdr=cat(3,Ni(:,:,1)*Theta(:,:,1)*Ni(:,:,5)'+2*Ni(:,:,3)*Theta(:,...
            :,1)*Ni(:,:,3)'-2*Ni(:,:,4)*Theta(:,:,2)*Ni(:,:,2)'+Ni(:,:,...
            5)*Theta(:,:,1)*Ni(:,:,1)'+2*Ni(:,:,2)*Theta(:,:,2)*Ni(:,:,...
            4)',Ni(:,:,5)*Theta(:,:,2)*Ni(:,:,1)'+2*Ni(:,:,3)*Theta(:,...
            :,2)*Ni(:,:,3)'-2*Ni(:,:,2)*Theta(:,:,1)*Ni(:,:,4)'+Ni(:,:,...
            1)*Theta(:,:,2)*Ni(:,:,5)'+2*Ni(:,:,4)*Theta(:,:,1)*Ni(:,:,...
            2)');
    case 'visc'
        gdr=cat(3,16*Ni(:,:,5)*Theta(:,:,1)*Ni(:,:,1)'+25*Ni(:,:,3)*...
            Theta(:,:,1)*Ni(:,:,3)'+9*Ni(:,:,1)*Theta(:,:,1)*Ni(:,:,5)'+...
            8*Ni(:,:,2)*Theta(:,:,2)*Ni(:,:,4)'+6*Ni(:,:,4)*Theta(:,:,2)...
            *Ni(:,:,2)',16*Ni(:,:,1)*Theta(:,:,2)*Ni(:,:,5)'+25*Ni(:,:,...
            3)*Theta(:,:,2)*Ni(:,:,3)'+9*Ni(:,:,5)*Theta(:,:,2)*Ni(:,:,...
            1)'+8*Ni(:,:,4)*Theta(:,:,1)*Ni(:,:,2)'+6*Ni(:,:,2)*Theta(:,...
            :,1)*Ni(:,:,4)')/9;
end

gdr=gdr(:);